Density matrix renormalization group algorithms with a single center site 
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We develop a correction to the density matrix used in density matrix renormalization group cal- 
culations to take into account the incompleteness of the environment block. The correction allows 
successful calculations using only a single site in the center of the system, rather than the standard 
two sites, improving typical computation times by a factor of two to four. In addition, in many 
cases where ordinary DMRG can get stuck in metastable configurations, the correction eliminates 
the sticking. We test the new method on the Heisenberg S — 1 chain. 



Since the density matrix renormalization group 
(DMRG) was devloped[HI3, it has gradually been ap- 
plied to more and more difficult systems, such as wide 
ladders and 2D clusters, and systems with long-range 
interactions. One of the problems arising in these sys- 
tems is the possibility that the simulation gets stuck far 
from the ground stateQ. Several approaches have been 
developed to overcome this problem, such as controlling 
the starting wavefunction through potentials or quantum 
numbers, with the controls later removed. Nevertheless, 
there has remained much room for improvement. In ID 
short-range systems, the standard DMRG finite system 
algorithm avoids convergence problems remarkably well 
because of the presence of the second center site in the 
block configuration. However, the extra site increases the 
computation time and memory requirements. An alter- 
native to utilizing the extra site, which works better in 
the more difficult cases, has not been available. In this 
paper we describe such an alternative method, which re- 
lies on a correction to the reduced density matrix in order 
to retain a broader variety of states. 

In the top panel of Fig. 1 we show the "superblock" 
configuration for the standard finite-system algorithm, 
where the lattice is divided into two large blocks, the 
system and the envirionment blocks, both with truncated 
bases, with two sites between them. The algorithm for a 
single DMRG step consists of finding the ground state for 
this "superblock" ; obtaining the density matrix for the 
system block plus site; diagonalizing this density matrix; 
and then changing basis to the most probable eigenvec- 
tors of the density matrix. This step replaces the system 
block, described by m states, by a block one site larger, 
but also described (approximately) by m states. One 
then shifts the dividing line between the system and en- 
vironment by one site, in order to add another site to 
the system block, and repeats the process. When the 
system block encompasses the whole system, the direc- 
tion is reversed and the roles of system and environment 
blocks are reversed. A sweep consists of one pass back 
and forth through the system. In a simple ID spin sys- 
tem one often obtains convergence to very high accuracy, 
e.g. an accuracy in the energy of order 10 -10 , with one 
or two sweeps through the lattice. 
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FIG. 1: Standard two-site DMRG method (top) and the 
single site method. 



In this description, it is apparent that one of the two 
sites in the center is crucial to the algorithm. The role of 
the other site is to increase the dimension and also the 
accuracy of the environment, particularly at the point 
where it connects to the system. One can leave out this 
extra site, i.e. use an environment block with the site 
already a part of it, as shown in the bottom panel of 
Fig. 1. This decreases the computation time for a step 
by roughly the number of states in a single site. How- 
ever, one finds that even in ID systems, the progress 
towards the ground state is much slower, and often stops 
altogether far from the true ground state. This can be 
understood in various ways. For example, suppose the 
ground state has total z component of spin S z = 0, and 
also suppose the environment block is poor and only has 
states with S z = 0. Then the renormalized system block 
will only have states with S z — 0, and no fluctuations 
in the spin will develop between the two blocks. In fact, 
any limitation on the quantum numbers present in the 
environment translates into a restriction on the states 
appearing in the renormalized system block. The distri- 
bution of states between various quantum numbers in the 
environment also translates directly to the renormalized 
system block. Note that if the environment block has m 
states, then the maximum number of nonzero eigenvalues 
of the density matrix is also m, and the number of states 
never increases unless states are added "artificially" de- 
spite having density matrix eigenvalues of zero. Simple 
fixes, such as adding extra random states with a larger 
range of quantum numbers, improve but do not fix the 
very poor convergence of the single site algorithm. 

The essential problem here arises when a particular 
fluctuation between the system and environment which 
should be present is not because the environment block 
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does not have the relevant states. Hence, the fluctua- 
tion is not represented in the density matrix and the new 
system block will not possess its relevant states for that 
fluctuation. Later, when the roles of system and envi- 
ronment are reversed, the relevant states again do not 
appear. In a ID system with short range interactions 
the extra environment site does a very good job of ensur- 
ing that the most relevant fluctuations are at least ap- 
proximately present the environment, so that subsequent 
sweeps can build in the fluctuations to high accuracy. 

In wide ladders or systems with longer range interac- 
tions, the addition of a single site to the environment is 
not always adequate. There may be missing fluctuations 
which are far from the extra site, and so are never built 
in. Even in these cases the extra site allows m to in- 
crease sensibly and as one lets m — > oo one obtains exact 
results. However, for practical values of m one may find 
unacceptably slow convergence. 

In this paper we describe an approximate correction 
to the density matrix to describe the key states which 
have been left out because the environment block is in- 
adequate. With this correction, the single site supcrblock 
configuration converges well. In addition, convergence in 
more difficult systems is dramatically improved, in either 
the single site or two site configurations. We present two 
different derivations of the correction, and give examples 
using the 5 = 1 Hciscnberg chain. 

We first give a simple, rough argument. Consider 
the power method for finding the ground state: iterate 
VVi+i = (1 — £ H)ip n , where £ is a small constant. As long 
as ip is no t orthogonal to the exact ground state, and 
e is small enough, the power method is guaranteed to 
converge to the ground state. Consequently, if the basis 
represents both ip and Hip exactly, and we minimize the 
energy within this basis, we expect exact convergence. 
The crucial point is the need to enlarge the basis to rep- 
resent Hip. Within the standard DMRG basis obtained 
from ip, after solving for the ground state, Hip = Eip, 
and nothing is changed by adding Hip to the basis. To 
go beyond the basis, we need to construct the parts of 
Hip as the basis is built up. The crucial terms of Hip 
come from the terms of H which connect the system and 
environment blocks. 

For the current supcrblock configuration, write the 
Hamiltonian in the form 

H = J2t a A a B a . (1) 

a 

Here the A a act only on the system block (including the 
site to be added to it), and the B a act only on the envi- 
ronment block (plus its site). All the terms which do not 
connect the blocks are contained in two terms of the sum 
which have either A or B equal to the identity operator, 
so that this form is completely general. (The other term 
in each case is the block Hamiltonian.) In order to put 
Hip into the basis, we need to target, in addition to ip, 



the terms A a ip for all a. Let the states of the system have 
indices s, p, and q, and the states of the environment e. 
The state A a ip can be written as 

£E^M*>l e >- (2) 

se p 

Targetting this wavefunction means adding into the den- 
sity matrix a term 

A/4, = a a A%iP pe r qe A a s , q * (3) 

epq 

where a a is an arbitrary constant determining how much 
weight to put into this additional state. The total con- 
tribution of all the terms is 

Ap = ]Ta Q i>i at (4) 

a 

where p is the density matrix determined in the usual 
way, only from ip. This is the form of the correction that 
we use, with a a — a ~ 1CP 3 — 1CP 4 . 

As a second derivation, we utilize perturbation theory. 
First, imagine that the environment block, but not the 
system block, is complete. We obtain the ground state 
exactly for this superblock, and then transform to the 
basis of density matrix eigenstates for the system block, 
and then also do the same for the environment block. 
Then the wavefunction can be written in the form 

\iP) = J2^s\L s )\R s ). (5) 

S 

The reduced density matrix is 

P = J2(R s \iP){ip\R s ) = J2 M 2 \Ls)(L.\ (6) 

s s 

Now consider the realistic case where the environment 
block is not complete. Assume the incompleteness takes 
the simple form that some of the \R S ) are missing, labeled 
s, whereas s arc present. Let P be a projection operator 
for the environment block P = ^2 S \s){s\, and take P — 
I— P. Let the unperturbed ground state, with energy Eq 
and density matrix p$, be obtained using the incomplete 
environment basis. We take as a perturbation the terms 
in the Hamiltonian which couple to the states s, namely 

H' = t a A a (PB a P + PB a P). (7) 

a 

The first order perturbative correction to the wavefunc- 
tion due to H' is 

\iP')=J2 t a (Eo-H )- 1 A a PB a \iP) (8) 

a 

where H = H-H'. 

In order to make progress we assume that each pertur- 
bation term A a PB a acting on the ground state creates a 
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set of nearly degenerate excited states, with average en- 
ergy E a . This assumption is equivalent to saying that the 
spectral function associated with each term is dominated 
by a narrow peak at E a . This significant approximation 
is reasonable because the correction to the density matrix 
is only used to enlarge the basis, to improve DMRG con- 
vergence. Correspondingly, we approximate (Eq — Ho)" 1 
as (Eq — Ea) -1 = l/e a . This gives 

W)^Y.^sY.- Aap z a \Ls)\R s ). (9) 

s a 

There are no first order corrections to the density matrix 
from \ip'), since P\tp) — 0. The lowest order correction 
to p can be written as 

Ap = J2 Ml E — — A a \L s )(L s ,\A ai M s , a , as (10) 

ss' aa' 

where 

M s , a , as = {R s ,\B a '^PB a \R s ) (11) 

Here if A is the unit operator, the term adds nothing 
to the basis. If B is the unit operator, M s > a / as vanishes. 
For the nontrivial pairs of operators A and B, this ma- 
trix element somewhat resembles a correlation function 
and it is natural to assume that the diagonal terms are 
dominant, where a — a' and s = s' . We expect the 
off diagonal terms a ^ a' to describe coherence between 
different perturbation terms which would tend to reduce 
the number of basis functions needed to describe the sys- 
tem block; therefore, ignoring the offdiagonal terms is a 
conservative assumption. Accordingly, we take 

M s i a r as ~ S ss 'S aa 'b a (12) 

This gives Eq. (4) with a a = b a \t a \ 2 /e^, and where we 
omit block-Hamiltonian terms. 

In practice, we take a a to be a small constant a inde- 
pendent of a. Construction of the correction to p take 
a calculation time for a single step proportional to m 3 
times the number of connecting terms, which is typically 
significantly smaller than the other parts of the DMRG 
calculation, although the scaling is the same. Larger val- 
ues of a introduce more "noise" into the basis, speeding 
convergence, but also limiting the final accuracy. Note 
that it is just as easy to apply the correction within the 
two-site method as the single-site method, which may be 
useful in some very difficult cases. We do not present 
results for this combination here. 

As a first test calculation, we consider the 5=1 
Heisenberg model 

H = J2SjS j+1 , (13) 

3 

where we have set the exchange coupling J to unity. The 
correction consist of the following: for each boundary site 
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FIG. 2: Error in the total energy for a 100 site Heisenberg 
spin-one chain, keeping m = 50 states per block, and using 
open boundaries terminated with 5 = 1/2 spins to remove 
the S = 1/2 end states (98 S = 1 sites + 2 S = 1/2's). 
The results are displayed for each half-sweep corresponding to 
reaching either the left or right end of the system. The two site 
method is the standard DMRG approach. The numerically 
exact energy was determined with the two-site method, using 
m = 200. 

1 of a block, i.e. a site directly connected to the other 
block, we add into the density matrix 

Ap = a(S+ P Sr + fir pS + + SfpSf). (14) 

For a chain with open boundaries, there is one site i; 
for periodic boundaries, there are two. One could argue 
that this expression should be adjusted with factors of 

2 between the z term and the other two terms, but this 
is not likely to make a significant difference. Note that 
the S + , S~ terms automatically increase the range of 
quantum numbers (i.e. total S z ) with nonzero density 
matrix eigenvalues. Figure 2 shows the convergence of 
the energy for a 100 site chain with open boundaries as 
a function of the sweep, keeping m = 50 states, relative 
to the numerically exact result obtained with m = 200 
and 10 sweeps. One can see the excellent convergence of 
the standard approach. The single-site method without 
corrections does not do too badly in this case, but still 
gets stuck significantly above the two-site energy. Adding 
the corrections, in this case with a = 10~ 4 , dramatically 
improves the convergence, making the single site method 
converge nearly as fast as the two site method. The two 
site method is roughly a factor of three slower than the 
single site method. Thus, even in this simple ID case 
where the standard approach works extremely well, there 
are advantages to using the corrected single site method. 

The results change significantly if we consider peri- 
odic boundary conditions. Here we consider the same 
superblock configuration as with open boundary condi- 
tions, but simply add in the connection to the Hamilto- 
nian between the first and last sites. There are better 
configurations for periodic boundaries, such as consider- 
ing it to be a ladder with the interchain couplings turned 
off except at the ends. These other configurations are 
superior only in the sense of improved convergence with 
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FIG. 3: Error in the total energy for a 100 site Heisenberg 
spin-one chain, with periodic boundaries. The number of 
states kept per block is indicated, and is the same for all 
three methods; four sweeps were made for each m. The cor- 
rection parameter a was taken to be 10~ 4 for sweeps 1-8, 
and 10 -6 for later sweeps. A somewhat slower convergence 
is visible for a — 10 -6 . The reference energy used was 100 
times the infinite energy per site, -1.401484038971(4) 0]. The 
corrected single site method using m = 4000 states gives a 
slightly lower total energy, due to exponentially small finite 
size effects, of -140.14840390392. 



the number of sweeps, not improved with respect to the 
number of states for a large number of sweeps. This naive 
configuration thus provides a difficult test for the single 
site method with corrections. In Fig. 3, we show the 
results for the same three cases as in Fig. 2. In this case, 
in the early sweeps, both uncorrected methods are stuck, 
ignoring the extra link between the first and last sites. 
The extra link eventually appears in the basis, but there 
is still sticking two or three times in higher energy states. 
In contrast, the corrected single site method never gets 
stuck and shows excellent convergence. 

A very useful DMRG technique is the extrapolation of 
the energy with the truncation error, i.e. the weight in 
the states which arc thrown out. If the truncation error 
were measured exactly, with a complete basis for the en- 
vironment, then the energy error would be proportional 
to the truncation error, allowing a linear extrapolation 
to zero truncation error. In practice, the apparent trun- 
cation error from the two site method may often be an 
underestimate, but one often finds that it is very con- 
sistent and still allows excellent extrapolation, even on 
fairly wide ladders. The truncation error within the cor- 
rected single-site method depends on a: as a — > 0, the 
apparent truncation error goes to zero and is unrelated 
to the exact truncation error. However, if a is not too 



small, linearity and excellent extrapolation are possible. 

Figure 4 shows results for the 100 site periodic system 
with a larger value of a, 10 -2 , suitable for extrapolation. 
The results show excellent linearity. The extrapolation 
gives -140.148416, off by 1.2 x 10~ 5 , whereas the sweep 
with to = 340 gave -140.148279, off by 1.2 x 10~ 4 . We 
have found that typically an order of magnitude improve- 
ment in the estimate for the energy is obtained by extrap- 
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FIG. 4: Error in the total energy for the system of Fig. 3 
versus the truncation error, with a — 10~ 2 . In this run two 
sweeps for each value of m, were made. The points shown are 
for m =80, 100, 120, 160, 200, 260, and 340. The line is a 
linear extrapolation, weighted with a standard deviation for 
each point assumed to be proportional to the truncation error 
at that point. 



olation in good cases; here we see similar improvement. 
In performing these extrapolations one always need to 
check the linearity for the system being studied. 

In summary, we have demonstrated a correction to 
the density matrix which allows the single-site DMRG 
method to converge well, and which improves the con- 
vergence dramatically for hard-to-converge systems. 
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